## This script produces Goodness of Fit diagnostics for "The Coevolution of Networks of Interstate Support, Interstate Threat and Civil War" 
## Since the estimation involves simulations, results may vary slightly from the published results.


library(reshape2)
library(network)
library(sna)
library(RSiena)
library(ggplot2)
library(gdata)
library(countrycode)
library(parallel)
library(ROCR)
library(PRROC)
library(miceadds)
library(sandwich)

rm(list=ls(all=TRUE))

#dir <- "[Directory Information Here]"
dir <- "C:/Users/rapiduser/Box/Relation creation/JOPreplication/"

setwd(dir)

(n.clus <- detectCores() - 1)

###############################################################################
##SienaGOF, Appendix Figure 10
###############################################################################
##Does not work with composition changes, so use only countries fully in data

#######Try with just a few years
keep(dir, n.clus, sure=TRUE)

nets.total <- read.csv("relation_network.csv", header=TRUE, row.names=1)
mons.total <- read.csv("relation_nodal.csv", header=TRUE, row.names=1)

## Set some parameters.

yr <- seq(2006, 2010)
yr.lo <- yr[1:(length(yr)-1)]

## List of names against which all data must be sorted. Also use this
## to build an empty matrix, to be used for merging data
id <- sort(unique(mons.total$ccode1))
n <- length(id)
emat <- matrix(NA, length(id), length(id), dimnames=list(id,id))


#########################################################
## Start with making the DV network matrices
#########################################################

##First the support matrices

for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","support_bin")]
  net <- acast(net, ccode1~ccode2, value.var="support_bin")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("S.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("S.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ss.NET", sep=""), lapply(objs,get))
names(Ss.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^S.NET."))

## Now the threat matrices
for (jj in yr) {
  if ( jj == yr[1] ) {
    sm <- 1
    objs <- vector()
  }
  net <- nets.total[nets.total$year==jj,c("ccode1","ccode2","threat_bin_3")]
  net <- acast(net, ccode1~ccode2, value.var="threat_bin_3")
  out <- emat
  out[rownames(out) %in% rownames(net), colnames(out) %in% colnames(net)] <- net
  ## Order rows/columns
  out <- out[order(id),order(id)]
  ## Rename
  assign(paste("T.NET.", jj, sep=""), out)
  rm(net,out)
  objs[sm] <- paste("T.NET.", jj, sep="")
  sm <- sm + 1
}; rm(jj)
assign(paste("Ts.NET", sep=""), lapply(objs,get))
names(Ts.NET) <- objs
rm(objs,sm)
rm(list = ls(pattern = "^T.NET."))


#########################################################
## Then do all the covariates for both networks
#########################################################

## Names of dyadic and monadic covariates, plus behavior variable
d.cons <- names(nets.total)[ !names(nets.total) %in% c("ccode1","ccode2","year","support_bin","threat_bin_3") ]
m.cons <- names(mons.total)[ !names(mons.total) %in% c("ccode1","year","intra_intens") ]
b.cons <- names(mons.total)[ names(mons.total) %in% c("intra_intens") ]

## Set the monadic covariates as NxK matrices, where K is number of years. Remember to
## subtract the last year for RSiena estimation.
for (k in 1:length(m.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr.lo) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(m.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr.lo
  out <- out[order(id),]
  assign(paste(m.cons[k],sep=""), out)
  rm(out)
}; rm(k)

#Now make behavior variable into matrix form, not substracting last year
for (k in 1:length(b.cons)) {
  out <- matrix(nrow=length(id),ncol=0)
  rownames(out) <- id
  for (i in yr) {
    tmp <- as.matrix(mons.total[mons.total$year==i,c(b.cons[k])])
    rownames(tmp) <- mons.total[mons.total$year==i,]$ccode1
    ## Bind on row names
    out <- cbind(out, tmp[match(rownames(out),rownames(tmp))])
    rm(tmp)
  }; rm(i)
  colnames(out) <- yr
  out <- out[order(id),]
  assign(paste(b.cons[k],sep=""), out)
  rm(out)
}; rm(k)

## Save dyadic covariates as lists of matrices
for (k in 1:length(d.cons)) {
  for (i in yr.lo) {
    cov <- nets.total[nets.total$year==i, c("ccode1","ccode2",d.cons[k])]
    ## Here's the adjacency matrix
    cov <- acast(cov, ccode1~ccode2, value.var=d.cons[k])
    out <- emat
    ## Add cov data to the full matrix with all actors
    out[rownames(out) %in% rownames(cov), colnames(out) %in% colnames(cov)] <- cov
    diag(out) <- 0
    ## Order rows/columns
    out <- out[order(id),order(id)]
    assign(paste(d.cons[k], i, sep=""), out)
    rm(out,cov)
  }; rm(i)
  out.list <- lapply(paste(d.cons[k], yr.lo, sep=""), get)
  names(out.list) <- paste(d.cons[k], yr.lo, sep="")
  assign(paste(d.cons[k], sep=""), out.list)
  rm(out.list, list=paste(d.cons[k], yr.lo, sep=""))
}; rm(k)

####################################
## Convert objects to RSiena objects
####################################

yrhi <- yr
yrlo <- yr.lo
nlo <- length(yrlo)
nhi <- length(yrhi)
N <- n

## Dependent networkS
assign("support.net", sienaNet(array(do.call("c", Ss.NET), dim=c(N,N,nhi))))
assign("threat.net", sienaNet(array(do.call("c", Ts.NET), dim=c(N,N,nhi)),allowOnly = FALSE))
assign("intra_intens.net", sienaNet(intra_intens, type="behavior"))

## Changing dyadic covariates, omitting the constant ones
dyad.effects <- d.cons[ !d.cons %in% c("contiguous") ]
for (x in dyad.effects) {
  assign(paste(x, ".net", sep=""), varDyadCovar(array(do.call("c", get(x)), dim=c(N,N,nlo))))
}; rm(x)

## Constant dyadic covariates
for (x in d.cons[ d.cons %in% c("contiguous")]) {
  assign(paste(x, ".net", sep=""), coDyadCovar(get(x)[[nlo]]))
}

## Monadic covariates
monad.effects <- m.cons
for (x in monad.effects) {
  assign(paste(x, ".net", sep=""), varCovar(get(x)))
}; rm(x)


## Create an RSiena network object

netdata <- sienaDataCreate(
  support.net,
  threat.net,
  intra_intens.net,
  
  mi1.net,
  energypc.net,
  mmp_change.net,
  dem.net,
  s_t_ipdistance.net,
  support_IGO.net,
  support_ipcloseness.net,
  support_mi.net,
  threat_mi.net,
  
  contiguous.net,
  idealpointdistance.net,
  IGOcommon.net,
  support_i_j_ipdistance.net
)

eff <- getEffects(netdata)

###Start with base model, and then add in information from support network

## Threat network equation
d.threat <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"),
  ".net", sep=""
)

m.threat_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)

m.threat_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.threat) { # Dyadic covariates
  eff <- includeEffects(eff, X, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_j) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, altX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.threat_i) { # Monadic covariates for threatening states
  eff <- includeEffects(eff, egoX, name="threat.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff <- includeEffects(eff, simX, name="threat.net", interaction1="dem.net", include=T, type="eval")

## Endogenous threat network effects
eff <- includeEffects(eff, recip, name="threat.net", include=T, type="eval")

#Now with civil war intensity as a dependent behavior variable
m.intra_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem"),
  ".net", sep=""
)

for (x in m.intra_i) { # Monadic covariates for threatened states
  eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1=x, include=T, type="eval")
}; rm(x)

## Add in information about the threat network for the behavior equation
eff <- includeEffects(eff, effFrom, name="intra_intens.net", interaction1="threat_mi.net", include=T, type="eval")


## Now Support network equation
eff1<-eff

d.support <- paste(
  c("contiguous",
    "idealpointdistance", "IGOcommon"
  ),
  ".net", sep=""
)
m.support_i <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)
m.support_j <- paste(
  c("mi1", "energypc", "mmp_change", "dem", "intra_intens"),
  ".net", sep=""
)


for (x in d.support) { # Dyadic effects
  eff1 <- includeEffects(eff1, X, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_j) { # Monadic effects for supported states
  eff1 <- includeEffects(eff1, altX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)

for (x in m.support_i) { # Monadic effects for supporters
  eff1 <- includeEffects(eff1, egoX, name="support.net", interaction1=x, include=T, type="eval")
}; rm(x)


## Do a political similarity measure
eff1 <- includeEffects(eff1, simX, name="support.net", interaction1="dem.net", include=T, type="eval")

## Endogenous support network effects
eff1 <- includeEffects(eff1, recip, name="support.net", include=T, type="eval")

## Include the cross-products
eff1 <- includeEffects(eff1, crprodRecip, name="support.net", interaction1="threat.net", include=T, type="eval")
eff1 <- includeEffects(eff1, crprodRecip, name="threat.net", interaction1="support.net", include=T, type="eval")

## Add in some higher order cross-network effects
eff1 <- includeEffects(eff1, inPopIntn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")
eff1 <- includeEffects(eff1, sharedIn, name="support.net",
                      interaction1="threat.net", include=T, type="eval")

## Add friend of my enemy effect
eff1 <- includeEffects(eff1, cl.XWX2, name="threat.net",
                      interaction1="support.net", include=T, type="eval")

##Now add information about the support groups
eff1 <- includeEffects(eff1, altX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff1 <- includeEffects(eff1, egoX, name="threat.net", interaction1="support_mi.net", include=T, type="eval")
eff1 <- includeEffects(eff1, effFrom, name="intra_intens.net", interaction1="support_mi.net", include=T, type="eval")
eff1 <- includeEffects(eff1, altX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff1 <- includeInteraction(eff1, altX, altX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff1 <- includeEffects(eff1, egoX, name="threat.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff1 <- includeInteraction(eff1, egoX, egoX, name="threat.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")
eff1 <- includeEffects(eff1, effFrom, name="intra_intens.net", interaction1="support_ipcloseness.net", include=T, type="eval")
eff1 <- includeInteraction(eff1, effFrom, effFrom, name="intra_intens.net", interaction1=c("support_mi.net", "support_ipcloseness.net"), include=T, type="eval")

#Estimate a model and use gofi function to extract goodness of fit
setwd(paste(dir, "/Output", sep=""))

tm <- Sys.time()
cl<-makeCluster(n.clus)
model <- sienaAlgorithmCreate(useStdInits=F, projname="Estimates.GOF.paper_1",
                              nsub=4, n3=5000, maxlike=F, firstg=.03,
                              modelType=c(support.net=1,threat.net=1))

output_1 <- siena07(model,data=netdata,effects=eff1,batch=TRUE,
                  verbose=FALSE,useCluster=TRUE,nbrNodes=n.clus,returnDeps=TRUE)

#plots, Appendix Figure 10
setwd(paste(dir, "/Figures", sep=""))

gofi_t_1 <- sienaGOF(output_1, IndegreeDistribution, verbose=TRUE, join=TRUE, 
                     varName="threat.net", cluster=cl)
summary(gofi_t_1)
pdf("GOF_ThreatIn-paper_1.pdf", width=5, height=3.5, paper="special")
plot(gofi_t_1)
dev.off()

gofo_t_1 <- sienaGOF(output_1, OutdegreeDistribution, verbose=TRUE, join=TRUE, 
                   varName="threat.net", cluster=cl)
summary(gofo_t_1)
pdf("GOF_ThreatOut-paper_1.pdf", width=5, height=3.5, paper="special")
plot(gofo_t_1)
dev.off()

gofb_1 <- sienaGOF(output_1, BehaviorDistribution, varName = "intra_intens.net")
summary(gofb_1)
pdf("GOF_Behavior-paper_1.pdf", width=5, height=3.5, paper="special")
plot(gofb_1)
dev.off()

goft_t_1 <- sienaGOF(output_1, TriadCensus, varName = "threat.net")
summary(goft_t_1)
pdf("GOF_TriadThreat-paper_1.pdf", width=5, height=3.5, paper="special")
plot(goft_t_1)
dev.off()




